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ABSTRACT 

For a magnetic flux tube, or indeed any flux, to emerge into the Solar corona from 
the convection zone it must pass through the partially ionised layers of the lower at- 
mosphere: the photosphere and the chromosphere. In such regions the ion-neutral 
collisions lead to an increased resistivity for currents flowing across magnetic field lines. 
This Cowling resistivity can exceed the Spitzer resistivity by orders of magnitude and 
in 2.5D simulations has been shown to be sufficient to remove all cross field current 
from emerging flux. Here we extend this modelling into 3D. Once again it is found 
that the Cowling resistivity removes perpendicular current. However the presence of 
3D structure prevents the simple comparison possible in 2.5D simulations. With a fully 
ionised atmosphere the flux emergence leads to an unphysically low temperature region 
in the overlying corona, lifting of chromospheric material and the subsequent onset of 
the Rayleigh- Taylor instability. Including neutrals removes the low temperature region, 
lifts less chromospheric matter and shows no signs of the Rayleigh- Taylor instability. 
Simulations of flux emergence therefore should include such a neutral layer in order to 
obtain the correct perpendicular current, remove the Rayleigh- Taylor instability and get 
the correct temperature profile. In situations when the temperature is not important, 
i.e. when no simulated spectral emission is required, a simple model for the neutral layer 
is demonstrated to adequately reproduce the results of fully consistent simulations. 

Subject headings: MHD - Sun: magnetic fields - Sun: chromosphere 



1. INTRODUCTION 

The emergence of new magnetic flux into the Solar corona is responsible for the formation 
of active regions. The accepted view is that the emergence of 0-shaped flux tubes through the 
photosphere is responsible for the formation of sunspots. It is clear therefore that the movement 
of magnetic flux from the convection zone up into the corona is one of the most significant drivers 
determining the structure of the corona. Furthermore the input of new flux into the preformed 
corona is often the trigger mechanism for dynamic coronal activity such as prominence eruptions. 
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flares and CMEs. Any attempt at a unified model of Solar activity must couple the magnetic field 
of the convection zone with that of the corona and hence a full understanding of the flux emergence 
process which connects them is essential. 

The problem with studying flux emergence is that it must couple sub-photospheric plasma with 
coronal plasma. In traversing this region the magnetic field moves from regions which are convec- 
tively unstable to convectively stable, through orders of magnitude changes in equilibrium density, 
and a rapid increase in temperature. The physics of each of these regions is therefore often dom- 
inated by different processes and analytical treatment of the whole emergence process is therefore 
limited. As a result this subject is now largely investigated by numerical simulations. A typical flux 
tube, assuming such a well defined structure e xists in the convection zone, must have sufficient twist 



to survive its transit of the convection zone ( Moreno-Insertis fc Emonet 



1996 



Porch &: Nordlund 



19981 ). It will then reach the photos phere where the buoyancy in s tability becomes act ive allowing 



the flux to escape into the corona (IMatsumoto &: Shibata 



1992 



Murrav et al.1 120061 ). Once the 



flux reaches the region of the chromosphere and corona, where the densit y drops by six orders of 



magnitude over a couple of megamet res in height, the flux tube expands (jMatsumoto et al 



1993 



Magara k. Longcopd l200ll : iFanI l200ll ). Below the photosphere the plasma beta j3 » 1 and the 
twist of the flux tube gives rise to a j x B force which is easily balanced by small changes in the 
much larger kinetic pressure terms. In the corona the plasma is characterised by /3 << 1 and any 
emerging j x B force cannot be balanced by kinetic pressure and thus the flux tube will expand 
rapidly into a configuration in which any residual j x B force is of the order of /?, i.e. the flux 
expands until the coronal field is force free. The initial twist in the emerging flux tube affects the 
emergence process and how close the coro nal field is to force free when it first reaches the corona 
(lAbbet fc Fisheill200.4 iMurray et al.ll2006l ) . This assumes that there is no overlying field with which 
the emerging flux can interact. Often such a field does exist and the inter action of this new fiux 
with existing magneti c structures can lead to complex, dynamic behaviour (lArchontis et al.l 12004 : 
Galsgaard et al.ll2005l ). The emerging magnetic flux, since it may have a non-zero j x B force, is 
also capable of lifting chromospheric material up into the corona. As the field expan ds, reducing 



the m agnetic forces, this heavier material may trigger the Rayleigh- Taylor instability (jlsobe et al 
2OO5I ). 



All of the works cited above have studied flux emergence by using the magnetohydrodynamic 
(MHD) equations appropriate for a fully ionised plasma. They make the further assumption that 
the parallel and perpendicular resistivities are equal whereas these differ by a factor ~ 2 for a fully 
ionised plasma. This factor of two is routinely ignored in MHD simulations as the resistivity used, 
for numerical reasons, exceeds the physical value by many orders of magnitude. For simulations 
of the coronal or convective regions of the Sun's atmosphere this is a perfectly valid approxima- 
tion. However, the photosphere and chromosphere are not fully ionised plasma due to their low 



( Cowling 


1957: 


Braginskii 


1965) 



resistivity, acts only on perpendicular currents, i.e. those flowing across the magnetic field, and 
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can be many orders of ma gnitude larger than the parallel Spitzer resistivity in the chromosphere 
([Khodachenko et al.l 12004 ). This dissipation of perpendicular curren ts by Cowling resistivity has 



been used to study the damping of MHD waves in the ch rom osphere (lGoodnia.nll2000l : iLeake et al 



20051 ) and flux emergence in 2.5D (jLeake &: Arbeii 120061 ) . In iLeake &: Arbeii (|2006l ) it was shown 



that the Cowling resistivity was sufficient to destroy all perpendicular currents in 2.5D flux emer- 
gence simulations. This is significant as it forces emerged flux to be force free as it traverses the 
chromosphere. The restricted geometry of 2.5D simulations prevents mass flow along the ignoreable 
direction and forces the emerg ence to be arcade like rather than bipolar. The aim of this paper is 
to return to the simulations in lLeake &: Arbeii (|2006l ) and study the effects of the partially ionised 
layers of the solar atmosphere on flux tube emergence in more realistic 3D geometry. 



2. EQUATIONS AND INITIAL CONDITIONS 
2.1. Equations 

The standard MHD equations, in Lagrangian form, are modified to include the effects of 
anisotropic current dissipation. For simplicity it is assumed that the atmosphere is composed 
entirely of hydrogen. The resulting set of equations apply to a single fluid and include the effects 
of partial ionisation through the neutral fraction ^„ = n„/ (rii -|-n„) where Un is the neutral number 
density and rii is the ion number density. 



£l = _ivP + ij AB + g + -V.S (2) 
Dt p p p 



(B.V)v-B(V.v) 

-VA(r?j||) - VA(7?J^) (3) 

tt: = V.v+-jj| +—j± 

Dt p p " p 

p T 

Where the parallel and perpendicular current vectors, jj| and respectively, are defined as 

- - w 

= — — 

and p is the mass density, P is the gas pressure, e is the internal specific energy density, v is the 
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centre of mass velocity of the fluid, B the magnetic field and g is gravitational acceleration. S is 



the stress tensor which has components Sij 



i<5,,V.v) 



and 



Ifdvi I dvj . 
2\'3x' "I" 'Bx'J 



Since the plasma is not fully ionised the total pressure P and specific internal energy e include 
the neutral fraction S,n through 



and 



+ (1 



nii 



(7) 



(8) 



^^m{l - 1) 

where ks is Boltzmann's constant, 7 is the ratio of specific heats, firn = fni/{2 — S,n) is the reduced 
mass and Xi is the ionisation energy of hydrogen. Equation H] can be used to numerically advance 
e but ^ „ is a function of temperature T so Equation [8] must be solved imp licitly for T which can 
then be used to specify P through Equation [71 For direct comparison with lLeake &: Arbeii (|2006l ) 
and all previous 3D flux emergence s imulations here we sim ply set /i^ = "^i /2 and ignore the Xi 
term in Equation (8). As discussed in lLeake &: Arbeii (|2006l ) this is unlikely to affect the emergence 
process through the chromosphere. The chromosphere is not in LTE and the radiation temperature 
and thermodynamic temperature cannot be assumed to be the same. The complete calculation of 
the ionisation state of hydrogen therefore requires the solution of the 3D radiative transfer and 
ionisation equations. To save time for the 3D emergen ce problem a simplified reduced model for 
^ri is used. This is based on a modified Saha equation (|Brownlll973l ) which can be solved for the 
steady state ionisation equatio n based on the local tem perature with the radiation temperature 
fixed at the photospheric value ([Thomas &: Athavlll96lh . 
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(9) 
(10) 
(11) 



where Tr is the temperature of the photospheric radiation field and w = 0.5 is its dilution factor. 
Using this equation, the ratio of the number density of neutrals to ions is given by 



m 



-1 + W 1 + 



and the neutral fraction = Pn/p is 



(12) 
(13) 



The dominant effect of the neutral atoms is how they modify Ohm's law and consequently lead 
to an an isotropic dissipation of current. The full derivation of the resistive terms in Equations (3) 
and (4) ( ICowlingj [19571 : [Braginskiil [1965[ ) shows that it is the collisions between ions and neutrals 
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which are most important. For a hydrogen plasma, and in the hmit when classical resistivity rj can 
be ignored, the Cowling resistivity is given by 

r]c = (14) 

with 



an = -?n(l-?n) \ ^m- (15) 

2 rUn V Trmj 

and Y^in = 5 x 10~^^m^ is the ion-neutral collisio n cross-section. When these formula are applied to 



model chromospheres (jKhodachenko et al.l 12004 ) it is found that the Cowling resist ivity r]r. can ex- 



ceed the classical resistivity r; by many orders of magnitude. From 2.5D simulations (jLeake &: Arber 



20061 ) it has been shown that the equations above include the dominant corrections to Ohm's law 
and that the neglect of the Hall term is valid for these flux emergence studies. Furthermore it was 
also shown that in the upper chromosphere the dominant contribution to Ohm's law is through the 
Cowling resistivity, not the advective term. 

The final term in Equation H] is designed to model all of the missing acoustic shock heating, 
radiative transport and thermal conduction terms. These terms would act to restore the equilibrium 
photosphere and chromosphere but are too computationally expensive, or simply unknown, and 
cannot be explicity included. As a first attempt at modelling the detailed heating and cooling 
terms omitted from Equation [4] we add a Newton cooling term — (e — eo)/T where r is the time- 
scale of the relaxation. The equilibrium specific energy density eo is chosen to be a function of 
the density p. The reasoning for this is related to the nature of these simulations. The buoyancy 
force drives magnetic field in the convection zone upwards into the photosphere, where the field 
then expands into the atmosphere above. Thus as a parcel of plasma from the convection zone 
of density p is moved upwards into the photosphere, its temperature should be relaxed to its own 
initial temperature, rather than the local plasma temperature, which is of different density. 



A form for the time-scale of this relaxation is required. For this the approach of lGudiksen &: Nordlund 



(|2005l ) is adopted. In simulating coronal heating they chose r to depend on some power of the den- 



sity ^ ^ 

r = 0.l(^) ■ (16) 
\PphJ 

so that at the relatively dense photosphere {p = pph) the time-scale is about 0.1s and is large 
enough that the effect becomes negligible in the sparse corona. 



2.2. Initial Conditions 

The modified MHD equations are normalised by division of the SI variables by photospheric 
values. The basic units are 



Lph = 150 km 
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Vph = 6.5km/s 
Pph = 2.7x lO-^kg/m^ 
Bph = 1200 G 

which gives the derived units 

tph = 23 s 
Tph = 6420 K 
Pph = 1.2x10'^ Pa 



From here on, unless stated, all quoted values are internal code variables and should be multi- 
plied by the above values to recover the SI variables. The differ ential equations (l )-(5) are advanced 



in time numerically using the Lagrangian remap code LareSd (jArber et al.ll200ll ). The physical do- 
main simulated extends vertically from -20 (3,000 km below the surface) to 130 (19,500 km above). 
The horizontal extend is 75 (11,250 km) about the centre of the domain, i.e. —75 > a; > 75 and 
—75 > y > 75. The z-axis is vertical, y across the tube and x aligned with the initial tube axis. 
Simulations have been run on 128^, 320^ and 512^ grids to check convergence. The computational 
grid was always uniform so that the minimum grid spacing used was Ax = 0.3, i.e. ~ 44 km. 

The anisotropy in the resistivity prevents the induction equation from being cast in simple dif- 
fusive form. In order to estimate the relative magnitudes of the resistive and the implicit numerical 
diffusion contributions to Equation (3) we consider the ID model equation dtB + vdxB = rjdxxB. 
For the second order accurate scheme employed here the leading order error term introduced in this 
model equation is of order vAx'^dxxxB so that if a typical scale length in the dynamic evolution 
is L this gives an effective numerical resitivity of vAx'^/L. In the upper chromosphere, where the 
Cowling resistivity is dominant, in nomalised units the maximum Alfven speed is 0.6. The worst 
case for numerical resolution corresponds to L ~ 1 as this would place three grid points across a 
slightly diffuse shock. This gives a normalised implicit numerical resistivity of about 0.06. Note 
that this estimate is based on the fastest phase speed and shortest gradient scalelength and thus 
a value of 0.06 represents an absolute upper limit on the numerical resistivity. The typical nor- 
malised r]c found in the simulations is of order 10, corresponding to a real magnetic diffusivity of 
9.75 X lO^m^/s, and is therefore larger than numerical resistivity. Note that this value of diffusivity 
corresponds to a magnetic Lundquist num ber of order 0.1. This is larger than that found for esti- 
mates based on the quite chromosphere in iLeake fc Arbe 3 jiooe). This is because chromospheric 



material expands as it is lifted due to flux emergence and the associated adiabatic cooling increases 
the Cowling resistivity. 

The initial stratification is a simple ID model of the temperature profile of the Sun, which 
includes the upper 3,000 km of the convection zone, the photosphere/chromosphere, the transition 
region, and the base of the corona. The temperature profile consists of a linear polytrope for the 
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convection zone with a vertical gradient at the critical adiabatic value 

dT _-f-lTdP 
dz 



(17) 



7 P dz 

The temperature in the photosphere and chromosphere is assumed to be constant at 1, as is the 
temperature in the corona at a temperature of 150. These two regions are connected by a transition 
region of width Wtr = 5. 

9 



T{z) 



Tph — 
Tph + 
tanh 



z, z < 



m + 1 

{tcor tph) 

2 

Wtr 



+ 1 



z > 



(18) 



(19) 



m 



-3Y is the adiabatic index for a polytrope, Zcor = 25 is the height of the corona, tph is the 



photospheric temperature and tc 



150. 



The density and pressure of the background atmosphere are found from solving the hydrostatic 
equation 

(20) 



dP_ 

dy 



-P9 



A magnetic tube is placed in the convection zone at z = — 10 with the profile 



Bo exp . 

qrB^ 



(21) 
(22) 



where r is the radial distance from the tube centre in the y, z plane. The strength of the field at 
the centre of the tube, Bq, is 5 and the the radius of the tube, a, is chosen to be 2. q is the amount 
of twist in the loop defined as 

Bx 



rB. 



(23) 



This is set to be the minimum require d to avoid fragmentation during the rise through the 
convection zone and is defined as \q\ = 1/a (jMoreno-Insertis &: EmonetlllQQq ). 



A choice must be made as to how to initialise the rise of the flux tube in the convection zone. 
It is thought that flux tubes forme d from the to roidal field in the tachocline remain connected to 
the large scale field by their roots (IZwannlll978l ). while the apex of the tube rises to the surface. 
As a result a flux tube which reaches the surface will be significantly 'bent' into an O-shape. In 
order to force the tube into this shape in these simulations, the centre is made buoyant while the 
ends are left in mechanical equilibrium. This is done by setting the pressure in the tube different 
to the field- free atmosphere {pq{z)) by pi(r) where 

dpi{r) 



dr 



JAB, 



(24) 
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so that the pressure gradient matches the Lorentz force. The density in the tube differs from the 
field-free density {po{z)) by pi{r) where 



a 



Pi 



X 
"A2 



(25) 



where a is used to scale the initial perturbation. In this paper, unless stated otherwise, a = 0.1. 
With this perturbation the centre of the tube, at x = 0, is buoyant whil e for x > A the tube is in 
mechanical equilibrium {pi = 0). The value of A is chosen to be 20, as in iFanI ([20011). 



2.3. Resistivity Models 

Three different models for the resistivity used in the flux emergence have been studied. The 
first was the fully ionised plasma model (labelled as FIP in later figures) in which ideal MHD was 
used. This is the same model used by all previous 3D flux emergence simulations and provides 
a benchmark against which the effects of Cowling resistivity can be measured. The second is the 
partially ionised plasma model (labelled as PIP in later figures) which solves the equations including 
partial ionisation, the Cowling resistivi t y and the Newton cooling as outlined above. This is the 



same model as used in iLeake Arberl (|2006l ). The final model is based on a simple model for 
partial ionisation effects in which a time independent perpendicular resistivity profile is fixed in a 
layer of the upper chromosphere. This model (labelled as Layer in later figures) has r] = and the 
Cowling resistivity fixed by 

r/c = 400^2 exp[-(z - 10) Vs] (26) 

in normalised units. This profile was chosen to closely match the resistivity observed in the chro- 
mospheric layers during simulations using the PIP model but still retains the magnetic field depen- 
dence. 



RESULTS 



The basic stages of th e flux emergence process are the same as in previous studies (e.g. iFan 



200 ll : lArchontis et al.ll2004l ) with the tube initially rising due to buoyancy. When the tube reaches 
the photosphere the buoyancy stops and the tube expands until sufficient flux has built up for the 
magnetic buoyancy instability to become important and the field expands through the chromo- 
sphere. The structure of the emerged field at t = 160 can be seen in Figure [H In previou s 2.5D 
studies of the effects of the partially ionised layers on flux emergence (jLeake &: Arberl |2006| ) it was 
possible to quantify these effects by calculating the integrated perpendicular and parallel currents 
as a function of height. These calculations were at the same time and only included flux inside the 
expanding envelope of flux. This simple measure of the effectiveness of the Cowling resistivity at 
removing perpendicular current is not possible in 3D due to the extra structure discussed below. 
Figure [2] shows the magnitude of j_L as a function of height along the line x = y = for all three 
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Fig. 1. — Field lines and flux for the partially ionised simulations at t = 160. The shaded contour 
plot, at the photopshere z = 0, shows vertical flux with positive flux shaded dark and negative flux 
shaded light. Dark field lines are those which connect to the initial equilibrium flux with r = 2.0 
while those shaded lighter correspond to fieldlines near the tube axis, i.e. r = 0.5. 



- 10 - 



resistivity models in the corona. While similar to the results in lLeake &: Arbeij (|2006l ) there are a 
number of significant differences. Firstly the results in Figure [2] all show a peak in at the top 
of the emerged flux. This is the expanding shock between flux free and emergi r ig flu x regions and 
was used to define the region inside which the integrated flux of lLeake &: Arbeii (|2006l ) was defined. 
More significantly the larger in the FIP model simulations cannot now be exclusively attributed 
to r]c removing this flux in the PIP simulations. The reason for this is that in 3D the chromospheric 
material raised into the corona in the FIP model is sufficient to trigger a Rayleigh- Taylor instability. 



Figure [3] shows the structure of the perpendicular current density in a vertical slice through 
the centre of the computational domain for the FIP and PIP models. The PIP model shows a slice 
through a symmetric, expanding shell of flux while for the FIP model there is a dropping central 
patch of enhanced j_L. This feature is due to the Rayleigh- Taylor instability, see discussion below, 
and as this leads to a bending and compression of fieldlines it also contributes to the net j^. Hence 
it is not possible to directly attribute the reduction of j_L in the PIP model, compared to the FIP 
model, shown in Figure [2] directly due to the dissipation of jj_ in the partially ionised chromosphere 
as some of the excess of j_L in the FIP model is created in situ in the corona by the Rayleigh- Taylor 
instability. 

The FIP model allows flux to emerge through the chromosphere with a larger j x B force 
than the PIP model in which the net jj_ is dissipated by Cowling resistivity. As a result the FIP 
model lifts more chromospheric material into the corona and is susceptible to the Rayleigh- Taylor 
instability. This is shown by the isosurface of density in Figure H] with the central density sheet 
dropping from the Rayleigh- Taylor instability clearly visible. Note that this density sheet is aligned 
with the magnetic field at the top of the emerging flux and hence, as can be seen from Figure [H 
is in a plane which crosses the tube axis. The x direction is therefore not an ignoreable coordinate 
for this Rayleigh- Taylor mode explaining its absence in the previous 2.5D work. 

A common feature of all flux emergence simulations using a fully ionised MHD model is that 
the rapid expansion of the emerging flux, once it reaches the low density corona, caused adiabatic 
cooling of the plasma to low temperatures. This can be seen for the FIP model in Figure [5] which 
compares the temperature as a function of height along a line at x = 10 and y = where the x 
coordinate is offset just enough to avoid the Rayleigh- Taylor instability induced density sheet. With 
the FIP model the temperature in the corona drops to 0.04 in normalised units, or ~ 250K. Including 
the Newton cooling/heating term with a relaxation time specified in Equation (16) only affects the 
temperature on the timescale of these simulations in a layer ~ 7 Mm above the photosphere and 
accounts for the ledge in temperature profile for both the Layer and PIP simulations between 
z = and z = 7. In the corona only the PIP simulations have a heating term due to the Cowling 
resistivity and this maintains the temperature there to about 0.7 times the photospheric value. 




Fig. 2. — Plot of the magnitude of j_L as a function of height along the line x = y = for all three 
resistivity models at t = 160. 
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Fig. 3. — Contour plots of the magnitude of j_L in a vertical slice through the computational domain 
in the x = plane. Results are for PIP model (top) and FIP model (bottom) at t = 160. 
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Fig. 4. — Volume rendering of density for the FIP model at t = 160. The volume contains values 
of density between 5 x 10~^ and 5 x 10~^. The color scaling on the cut through the volume shows 
higher density as darker shades. The domain has been split along the axis of the initial equilibrium 
flux tube to show the density structure inside the expanding shell being lifted by the flux emergence. 
The dominant Rayleigh-Taylor spike can be seen dropping at around a; = in the y = plane. 
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Fig. 5. — Temperature as a function of height for ah three resistivity models at t = 160. 
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CONCLUSIONS 



Previous 2.5D simulations ( Leake &: Arb had shown that including the partially ionised 



layers of the Solar atmosphere changed the dynamics of magnetic flux emergence by introducing a 
Cowling resistivity. This dissipated perpendicular current density and lead to a force free coronal 
field (inside the expanding region of emerged flux but ignoring the perpendicular current density 
at the interface of emerged flux and field free corona). This simple picture cannot be supported in 
such a clear way for the full 3D simulations presented here. In 3D the most pronounced difference 
between partially ionised and fully ionised simulations is that the fully ionised simulations raise 
more chromospheric material. For the initial conditions used in this paper this meant that the fully 
ionised model became unstable to the Rayleigh- Taylor instability while the partially ionised model 
remained symmetric with no signs of instability. The presence of the Rayleigh- Taylor induced 
perpendicular current density means that it is not possible to assess the affect of the Cowling 
resistivity on the amount of perpendicular current density emerging into the corona in isolation. 

The primary result of these simulations is therefore that the inclusion of the Cowling resistivity 
affects the amount of chromospheric material uplifted into the corona. The process of removing 
jj_ allows most of the field to move through the chomosph eric plasma rathe r than lifting it, thus 



reducing the amount of mass uplifted. Recent publications (jlsobe et al.ll2005l ) have suggested that 
the onset of the Rayleigh- Taylor instability during flux emergence may be the cause of coronal 
loops. This result may depend critically on the neutral hydrogen but this was absent from all 
previous flux emergence simulations. 

Often in coronal physics one is only concerned with the emergence of magnetic flux and its 
structure. In such situations this paper has shown that a greatly simplified model of the Cowling 
resistivity is capable of reproducing the results of the full PIP simulations, except for the tempera- 
ture profile. Since the present, and indeed all previous, flux emergence simulations have omitted a 
full treatment of thermal conduction, radiation effects and coronal heating it is unlikely that such 
simulations can achieved accurate temprature estimates for emergence. A practical conclusion from 
this work is therefore that the minimum physics required to obtain credible fleld structures in flux 
emergence is the Layer model presented in Equation [26l This has the advantage of being easy to 
include into any code and yet accurately predicts the correct magnetic field structure and uplifting 
of chromospheric material. 



This work was funded in part by the Particle Physics and Astronomy Research Council. The 
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and Warwick University's Centre for Scientific Computing. 
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